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From three-dimensional linearized hydrodynamic equations, it is found that the heat 
conductivity is proportional to [L^ / {L^LDY^^ , where Lx, Ly and are the lengths of the 
system along the x, y and z directions, and we consider the case in which ^ Ly^L^. The 

necessary condition for such a size dependence is derived as = L^j ['n}^'^ iJ'J'^ iJ'J^') ^ 1, 
where </> is the critical condition parameter and n is the number density. This size dependence 
of the heat conductivity has been confirmed by molecular dynamics simulation. 



§1. Introduction 

Fourier's law of heat conduction is a relation between the heat flux J and the 
temperature gradient. It is given by J = — AVT, where A and T are the heat 
conductivity and the temperature, respectively. Although it is widely believed that 
Fourier's law is realized in many situations, recent numerical results suggest that the 
heat conductivity diverges in low-dimensional systems in the thermodynamic limit, 
while it is convergent in three-dimensional (3D) systems.^)' 

Narayan and Ramaswamy^^ have found that the heat conductivity is propor- 
tional to L^/^, with the system size L. They derived this result from hydrodynamic 
equations with thermal fluctuations in the ID limit. This size dependence of the heat 
conductivity is also applicable to ID chains. These results have been confirmed in 
numerical simulations of the Fermi-Pasta-Ulam (FPU) chain, a ID gas modeF^'^-* 
and single-walled carbon nanotubes (SWNT).^)' ^'^^ 

In real experiments, of course, it is difficult to realize actual one- or two-dimensional 
systems. However, it is easy to prepare quasi-one-dimensional (QID) systems in 
which the role of one direction is dominant, and the length of the system along that 
direction is much greater than those along the other directions. Similarly, we can 
construct quasi-two-dimensional (Q2D) systems, in which two directions are domi- 
nant, the lengths of the system along those directions are much greater than those 
along the other direction. In this paper, we demonstrate that heat conductivities in 
QID systems diverge as A ~ (La;/(LyL^))^/^, where is the length of the system 
along the x-direction, and we have 3> Ly,Lz. The necessary condition for this 
anomalous behavior of the heat conductivity is derived as 

where n is the number density, is the number of particles in the system, and 
(j) is the critical condition parameter. We also find that Q2D systems diverge as 
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X ~ y(ln(L^ + Ly))/Lz, where and Ly are much larger than L^- The necessary 

condition to obtain this behavior is ((L^ + Ly) / (riLl)) -|- L^/(toCo)] 1, 

where to is the mean free time and cq is the velocity of sound. 

§2. Derivation of the long-time tail from the hydrodynamic equations 

The divergence of the transport coefficient in the thermodynamic limit originates 
from the long-time tail of the time correlation function, which has been confirmed 
by molecular dynamics (MD) simulations with hard spheres. For the theoretical 
calculation of the long-time tail near the equilibrium state, the mode-coupling the- 
ory^^) is a powerful tool. Amongst several methods for calculating the long-time tail 
in the mode-coupling theory, here, we adopt the hydrodynamic approach developed 
by Ernst et al}^^'^"^^ 

From the Green-Kubo formula, ' the heat conductivity A can be calculated 

as 

X = (3T-^ lim f'dtCxit), (2-1) 

where C\{t) is the time correlation function of the heat flux J and (3 = (kBT)^^ 
with the Boltzmann constant ks- Here, we define the x-component of the heat flux 
J as 




where Vix is the x component of the velocity of the particle z in a fluid of N particles, 
m is the mass of a particle, rij^x is the x component of the displacement vector 
between particles i and j, and (^(r) is the intermolecular pair potential. From this 
point, we consider only the kinetic parts, Jk, of J, given by 

J^^E (^-" - I'^bt] V.X, (2-3) 

i ^ ^ 

and the corresponding correlation function. Through simulations, as reported in ^ 
we have verified that the error introduced by ignoring the remaining terms in (j2-2p 
is negligibly small. 

Although the Green-Kubo formula can only be proven for infinite systems, if we 
apply it to finite systems, we need to introduce upper bound of the integral tc, which 
represents the typical transit time, tc ~ L/cq, without taking the limit — > oo.^^ 

To calculate the time correlation function, we need to solve the d-dimensional 
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linearized hydro dynamic equations 
dn{r, t) 



dt 
du{r,t) 



- — nV ■ u{r, t), 

=i/V^M(r, t) + {Di - z^) V (V • u{r, t)) 



dt 

cl fVn{r,t) A (2-4) 

° ' ^ ' ^ +apVT{r,t) > 



7 V ^ 
dT{r,t) 7-1 



V • ii (r, t) + 7L>TV^r (r, . 



9i aj; 
Here, u is the velocity field and 



7 = — , Co 



m \ dn J rp 



A 2 (d - 1) r? + 

= — , Ut = — , Ui = , 

nm ncp anm 



(2-5) 



where p is the pressure, and Cp and cy are the heat capacities per particle at constant 
pressure and at constant volume, respectively. Also, ij and ^ are the shear viscosity 
and bulk viscosity, respectively. 

To solve the set of equations ()2-4p . Ernst et al.^^^ use the Fourier transform 
with respect to space. In this paper, we use the Fourier series of the hydrodynamic 
variables with respect to space, given by 



nq(t) = [ dr ,e-''^-'^6n{r,t) , 
Jv 



V 

(2-6) 



where 5n = n{r, t) — n is the deviation from the equilibrium value, and the system is 
a rectangular solid {L^,, Ly, Lz}, qi = 2'Kki/Li {i = x, y, z), with integer ki. Similarly, 
we introduce the Fourier components Tq{t) and Uq{t), which are the components of 
6T (r, t) = T {r,t) — T and u{r,t), respectively. 

In d-dimensional Fourier space, Eq. (j2-4|) is converted into 



dt 

dUq 

dt 



mq • Uq, 

uq^Uq - {Di -v)q{q- Uq) 



J \ n '■I 



— T^- = iq ■ Uq - jDrq Tq. 

ot an 



These equations can be reduced to a diagonalizable (d+2)-dimensional matrix equa- 
tion. Here, we introduce the "hydrodynamic modes", defined as the eigenfunctions 
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of the form exp{—ujqt). The quantities ojq in the small q (= |q|) limit are calculated 
as the d — 1 viscous modes [tOq = vq^ {j = " i d — 1)], the two sound modes 
{uj^ = ±icoq +lrs), with 

r, = A + (7 - 1) Dt, (2-8) 

and the heat mode = D^q^)- The Fourier series of the hydrodynamic variables 
can be expressed as linear combinations of the eigenvectors of the matrix equation. 

After setting the initial values used by Ernst et al.,^^^ C\{t) can be calculated. 
When the system is sufficiently large, the most dominant part of the long-time tail 
of Cx{t) is the viscous-heat (VH) term, given by 



(2-9) 



with qi = qi/q {i = x, y, z) and the sound-sound (SS) term, 



(2-10) 



Note that Cx{t) decays as f-'i/2.i3),i4) 

In a QID system, however, the SS term is dominant, as seen from the following. 
First, let us assume that Ly ~ and that is sufficiently larger than Ly and L^. 
The exponential term in the SS term in Eq. (|2-10|) is 



exp 



-47r2A 



kl kl 



r2 ^ r2 ^ r2 



(2-11) 



LI + LI + Ll^L, 



When time is sufficiently large, i.e., t ^ tc ^ L/cq with L = 

it is found that the {ky,kz} = {0,0} term is dominant in Eq. (|2-10|) . Because we 
consider QID systems, we assume that the relation given by 



:L^ > 1 



L2 cqLI' 

is satisfied. Therefore Eq. (j2-10p can be rewritten as 



1 1 

LyLz Lx 



E 



exp 



A-2 ■ 



(2-12) 



(2-13) 



If Fst/Ll. ~ Fs/^cqLx) is sufhciently small, this summation can be reduced to the 
integral form 



1 f"^ 
LyLz Jo 



AnLyLz 



vr 

TTt' 



(2-14) 
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which imphes that C^^{t) decays as 

By contrast, the {ky,kz} = {0,0} term in the VH term of Eq. (j2-9p cannot be 
dominant. Instead, the {ky,kz} = {±1,0} and {0,ibl} terms are dominant, and 
Eq. (j2-9p can be written in the integral form 



LyLz JQ 



dx 



1 + L2x2 



exp 



-A-n'^{v + Dt)- 



1 + Lix 



LI 



< 



LyLz 



exp 



1/2 



/>oo 

/ dx exp [-47r2 {v + Dt) x^t] , (2- 

JO 



15) 



vr 



ATvLyL^y {u + DT)t 



which reveals that C^^{t) decays exponentially with time. Thus, C^^{t) decays 
much faster than C^'^{t). From Eq. (j2Tp . if is seen that the heat conductivity can 
be calculated as 



A ~ PT-^ / ' dt {Cf (t) + C7]f^ (t)) ~ PT-^ / ' dtC^^it). 



(2-16) 



From Eqs. (j2-14p and (j2-16p . we need to estimate the system size dependence of 
Fg in order to calculate A. Therefore, we need to calculate the long-time tail of the 
shear viscosity in QID systems. The time correlation function of the xy component 
of the shear stress and the corresponding viscosity are calculated as 



ri^y = ^f^^dtc^^y{t), 
ciy{t)^{r,y{t)r,ym, 



•Xt)Viy{t) — - / ^^tj.a; 



(2-17) 
(2-18) 

(2-19) 
(2-20) 



We consider only the kinetic parts, J^'^^ , of J^^ and the corresponding correlation 
function. From the similar calculation of Cx{t), if is found that the dominant term 
of C^{t) for large is the viscous- viscous (VV) term from the yz component of the 
shear stress, given by 

^^(t) ~ + 24^41] eM-^Wt] ~ L.y^. (2-21) 



Because the kinetic viscosity is proportional to r] in Eq. ()2-5p . we derive r/ from 
Eq. (l2Tr|) as 
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We assume that the bulk viscosity ^ satisfies the relation 

which is realized in most normal fluids. Therefore, Fg can be expressed as 

/ r \i/3 

r,r^CiX + C2^j^J , (2-24) 

where Ci and C2 are constants that are independent of time and the system size. 
Thus, from Eqs. (j2-14p . ()2-16p and (j2-24p . we can derive the following self-consistent 
equation for A: 

CiX' + C2(^Y\'-^. (2-25) 

This equation indicates that the system size dependence of A is given by 

/ r \i/3 

The necessary condition for Eq. ()2-26p to be valid is derived from Eq. (|2-12p as 

(2-27) 

The heat conductivity and shear viscosity in Q2D systems can be derived cor- 
respondingly as 

— ' ^^'^^^ 

where <C Lx,Ly. The necessary condition for this anomalous behavior is derived 
as 



y In ^ » 1. (2-29) 

nL^ toCo 

S3. MD simulation 



In order to check the validity of Eq. (j2-26p . we have performed 3D MD simula- 
tions with hard spheres in QID systems. ^^•'"^^^ It is noted that a logarithmic system 
size dependence of the heat conductivity in 2D systems has been confirmed with MD 
simulations with hard disks. ^^-^ It should also be noted, however, that it is difficult 
to distinguish logarithmic behavior from power law behavior in 2D simulations. 

In our simulations, the system is confined in a box of size >i Ly x L^, as shown 
in Fig. [H The walls perpendicular to the y and z directions are perfectly refiecting 
walls and the walls vertical to the x direction are perfectly thermalizing walls at 
temperatures Ti and T2, respectively. At the 'thermal' wall of temperature T^, a 
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Fig. 1. The system studied in the simulation. The walls perpendicular to the y and z directions 
are perfectly reflecting and those perpendicular to the x direction are perfectly thermalizing. 
These 'thermal' walls are at Ti and T2, respectively. 



particle is reflected with a new velocity v at random. The probability distribution 
function for v is given by 



1'k\Vx\ 



{2TTkBTu, 



■ exp 



V 



2k bT^ 



(3-1) 



The MD data were taken for the systems with Ly = = 0.1 — 3.2 and L^/Ly = 
2 — 16384, where the unit of length is the diameter of a particle and Li = Li — \ 
is the range within which the center of the particle can move. For the system with 
Ly = 3.2, data were taken only for L^/Ly = 2 — 4096. The packing fraction was 
fixed approximately to 0.055. The ratio of the temperature difference between two 
thermal walls was T1/T2 = 2. 

As the initial state, the molecules were arranged so as to realize a constant 
pressure under the linear profile of T; i.e. nT is initially uniform in the system. The 
initial velocities of the particles were chosen from the Maxwellian distribution with 
T varying linearly from Ti to T2. To obtain a stationary initial state, we simulated 
6 X 10^ collisions per particle which started from the initial state. We carried out 
the simulation until 4 x 10^ collisions per particle were realized from the stationary 
initial state. The heat flux J was calculated using Eq. (|2-2p . We found that in the 
simulation, the difference between J and Jk is less than 0.1%. The heat conductivity 
is calculated as A = —JLx/^T, where AT is the difference between the temperatures 
near the thermal walls at Ti and T2, respectively, and Lx is the distance between 
the points at which the temperatures are given. 

The result of the simulation is shown in Fig. [2j There, it is seen that the heat 
conductivity seems to obey Eq. (|2-26|) for Ly > 1.8, while such behavior is not 
observed for smaller Ly. We believe that this difference is caused by the limitation 
on the collision angle in the case that Ly is much smaller than the diameter. For 
small Ly, as in Fig. [3l most of the collisions of the particles take place with v'l^ ~ V2x 
and V2x — vix , and the changes of the velocities in the other directions are small in 
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Fig. 2. The size dependence of the heat conductivity in the QID system. The x-axis is L^/Ly, 
with the unit of length being the diameter of the hard spheres. The circles, crosses and open 
squares indicate the heat conductivities, where Ly is 1.8, 2.6 and 4.2, respectively. The solid 
line represents a function proportional to x^^^. 



Fig. 3. A collision of particle 1 and particle 2 for small Ly. Here, vi and V2 are the velocities of 
particles 1 and 2 before the collision, respectively, and v'l and v'2 are the velocities of particles 
1 and 2 after the collision, respectively. 



each collision. Thus, most of the colhsions cause only exchanges of the velocities of 
the particles in the x-direction, and the heat transfer behaves like that in the case 
of ballistic transport. In this situation, the linearized hydrodynamic equations (j2-4p 
for small Ly cannot be used. Contrastingly, for the system with Ly = 1.8 and 
Lx/ (LyLl) > 10.0, the most probable value of the exponent a with is 0.30 ±0.01, 
which is close to 1/3. 

Although the results of our theory and simulation are consistent, we need to 
simulate larger systems to confirm the relation (j2-26p . There is an inconsistency in 
the analytical and numerical treatments. Indeed, because the linearized hydrody- 
namic equations (j2-4p were adopted in the case that the system size is larger than 
the mean free path, but this conditions is not satisfied in our simulation, where the 
widths Ly and L^ are smaller than the mean free path by an amount on the order 
of 7 diameters. 




X 
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Fig. 4. The size dependence of the heat conductivity in the QID system. The s-axis is Lx, with the 
unit of length being the diameter of the hard spheres. The circles, crosses and open squares indi- 
cate the heat conductivities, where Ly is 1.8,2.6 and 4.2, respectively. The solid line represents 
a function proportional to x^^^. 




§4. Discussion 

We should note that Shimada et al?^ and Ogushi et al.^^ also calculated the heat 
conductivity using 3D MD simulations in QID systems consisting of hard spheres 
and Lennard-Jones molecules, respectively. Their simulations suggest that the heat 
conductivity does not diverge in the limit ^ oo. However, the system sizes in 
their simulations do not satisfy the necessary condition 0^1. In their simulations, 
(j) is 0.03 — 0.3 with Ly = 20.0 and 0.3 — 2 with Ly = 4.0, respectively, where the unit 
of the length is the particle diameter. By contrast, our simulation is more extensive 
in the system cj) is 1.8 — 2 x 10*^ for Ly = 1.8. This difference might be the cause of 
the qualitative difference in the behavior of the heat conductivity. 

In numerical simulations of SWNT, which can be regarded as a QID system, 
A ~ has been observed. In Murayama's paper, it is reported that the 
heat conductivity behaves as A ~ Lf^"^, with the tube length Li for (5,5) SWNT. 
However, for a wider tube diameter, (10, 10) SWNT, the increase of A is suppressed. 
It might be possible to understand this result in terms of the necessary condition 
for the anomalous behavior of the heat conductivity, 0^1. This behavior is also 
consistent with that found in our simulation. As we can see in Fig. [H anomalous 
behavior of the heat conductivity, a ~ 1/3, with A ~ L", is observed for L^ > 10^, 
with Ly = 1.8 where (p is larger than the critical value, 4>c = 7.1 x 10^. By contrast, 
for a wider tube diameter, Ly = 4.2, the increase of A is suppressed. From the 
necessary condition (j2-27p . the anomalous behavior q ~ 1/3 may be observed for 
Lx > 8.4 X 10'^ for the wider tube diameter Ly = 4.2, for which cj) > cpc satisfied. 
Similarly, for (10, 10) SWNT, anomalous behavior of the heat conductivity might be 
observed for a longer tube, because (/> » 1 is satisfied in that case. 

Recently, Shiba et al?^^ observed that the time correlation function of the heat 
flux is proportional to i-0-98±0-25^ g^j^j i^j^g heat conductivity exhibits logarithmic 
divergence with the system size in three-dimensional FPU-/3 lattice systems. They 
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also measured the size dependence of the heat conductivity in QID FPU-/3 lattices. 
Figure 4 in their paper suggests that when the system is sufficiently close to a QID 
system, the heat conductivity is proportional to L^/^. 

In conclusion, we derived relation according to which the heat conductivity 
is proportional to (Lx/LyLl)^/^ in QID systems and (In (L^ + L^))/Ly in Q2D 
systems, and the necessary conditions L'^. / {n^ LyL^) » 1 and 

((L^ + L^)/(nL^)) ln[y^Z2~+~L2/(tgCo)] 1, respectively, for the anomalous behav- 
ior of the heat conductivity from linearized hydrodynamic equations. This behavior 
in QID systems has been confirmed through comparison with MD simulations. 
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